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Abstract 

Previously, we have presented a methodology to extend canonical 
Monte Carlo methods inspired on a suitable extension of the canonical 
fluctuation relation C = /3 2 compatible with negative heat capac- 

ities C < 0. Now, we improve this methodology by introducing a better 
treatment of finite size effects affecting the precision of a direct determi- 
nation of the microcanonical caloric curve fi{E) = dS (E) /dE, as well 
as a better implementation of MC schemes. We shall show that despite 
the modifications considered, the extended canonical MC methods possi- 
bility an impressive overcome of the so-called super- critical slowing down 
observed close to the region of a temperature driven first-order phase tran- 
sition. In this case, the dependence of the decorrelation time r with the 
system size N is reduced from an exponential growth to a weak power- 
law behavior t(N) tx N a , which is shown in the particular case of the 2D 
seven-state Potts model where the exponent a = 0.14 — 0.18. 



1 Introduction 



In a previous paper [I], we have proposed a methodology that enables Monte 
Carlo (MC) methods based on the Gibbs canonical ensemble: 

dp c {E\p B ) = -± n -exp{-0 B E)Sl(E)dE (1) 
* (Pb) 

to account for the existence of an anomalous region with negative heat capacities 
C < [2j [3l [H [6j [7] and to overcome the so-called supercritical slowing down 
observed near of first-order phase transition [5]. This development is inspired 
on the consideration of a recently obtained fluctuation relation [§J [TUJ [TT] : 

C = 1 (6E 2 ) + C (5j3 u 5E) , (2) 

which appears as a suitable extension of the known canonical identity C = 
1 (SE 2 ^ involving the heat capacity C and the energy fluctuations (5i? 2 ). This 
last expression accounts for the realistic perturbation provoked on the internal 
state of a certain environment as a consequence of the underlying thermody- 
namic interaction with the system under study, which is characterized here in 
terms of the correlation function (Sp^SE) between the system internal energy E 
and the environment inverse temperature (3 U . While the canonical ensemble (QJ 
dismisses the existence of this feedback effect (due to the constancy of the in- 
verse temperature /3b associated with this ensemble) and it is only compatible 
with positive heat capacities, the general case possibilities to access anoma- 
lous macrostates with negative heat capacities C < as long as the condition 
(5P U 5E) > 1 is obeyed. 

The incidence of non- vanishing correlated fluctuations (6f3 u 5E) can be easily 
implemented in MC simulations. Roughly speaking, the extension of canonical 
MC methods is achieved by replacing the canonical inverse temperature j3s with 
a variable inverse temperature j3 u (E). The resulting framework constitutes 
a suitable extension of the Gerling and Hiiller methodology on the basis of 
the so-called dynamic ensemble [12] , where the microcanonical curve /3 (E) = 
dS (E) I dE and the heat capacity C (E) can be estimated in terms of the energy 
and temperature expectation values (E) and (ft J) as well as their fluctuating 
behavior described by Eq.@. This method successfully reduces the exponential 
divergence of decorrelation time r cx exp (XN) with the increase of the system 
size N of canonical MC methods to a weak power-law divergence r cx N a , 
with a typical exponent a ~ 0.2 for the case of 2D ten-state Potts model 
[1 . By combining this type of argument with cluster algorithms, one obtains 
very efficient MC schemes that constitute attractive alternatives of the known 
multicanonical method and its variants [5]. In this work, we shall improve 
the present methodology to consider the existence of finite size effects that 
reduces the precision of a direct determination of the microcanonical caloric 
curve /? (E) — dS (E) /dE, as well as a better implementation of MC schemes. 
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2 Methodology 



2.1 Overview 

The simplest way to implement the existence of non-vanishing correlations 
{5f3 u 8E) corresponds to a linear coupling of the environment inverse tempera- 
ture with the thermal fluctuations of the system energy: 

/3 U (E) =(3 e + X5E/N, (3) 

where A is a coupling constant that appears here as an additional control pa- 
rameter. While the case with A = corresponds to the canonical ensemble ([1]), 
where /3 U — const, in general, the constancy of the inverse temperature /3 e can 
be only ensured in the average sense, j3 e — (/3J). Eq.Q can be substituted into 
Eq.Q to obtain the following results: 

(Ai?)2 = wwcr-y {A ^ 2 = ^^WTa' (4) 



where Ax = \J (5x 2 ) denotes the thermal dispersion of a given quantity x. 
Since A/3 W and AE should be nonnegative, one arrives at the following stability 
condition: 

f3 2 N/C + \>0. (5) 

For A = 0, one derives the ordinary constraint C > that emphasizes the 
unstable character of macrostates with negative heat capacities C < within the 
canonical description. However, these anomalous macrostates can be observed 
in a stable way in a general situation with A =/= when this control parameter 
satisfies Eq.©. By assuming an extensive character of the heat capacity C ~ N 
in short-range interacting systems, the energy dispersion AE grows with the 
increase of the system size N as AE ~ >/~N, so that, the dispersion of the 
energy per particle e = E/N behaves as Ae ~ l/y/~N. Since 5/3 u = XSE/N in 
the case of the ansatz <j3j> , the dispersion of the inverse temperature also behaves 
as A/? w ~ 1/y/N. Thus, the present equilibrium situation constitutes a physical 
scenario that ensures the stability of macrostates with negative heat capacities 
C < with the incidence of small thermal fluctuations. 

This type of equilibrium situation is schematically represented in FIGfTJ 
Here, it is shown the typical microcanonical caloric curve (3 (E) = dS (E) /BE 
corresponding to a finite short-range interacting system that undergoes a first- 
order phase transition, as well as the energy distribution function p (E) associ- 
ated with the thermal coupling of this system with an environment with inverse 
temperature f} u (E). The intersection points E e derived from the condition of 
thermal equilibrium: 

0{E a ) = u (E e ) (6) 

determine the position of maxima and minima of the distribution function p (E). 
Note that the linear ansatz (|3]) can be considered as the first-order approxima- 



3 





C<0 

i> 


/ 
P„( E > 






' — • / e 








P / 










P(E) 

/ 


P(E) 

\ 



Figure 1: Schematic representation of the typical microcanonical caloric curve 
/3 (E) = dS (E) j dE associated with a finite short-range interacting system 
undergoing a first-order phase transition, as well as the energy distribution 
function p (E) resulting from the thermal coupling of this system with a certain 
environment with inverse temperature /3 U (E). For further explanations, see the 
text. 



tion of the power expansion of a general dependence /3 U (E) : 



oo 



ME) =/3 e + 2^a n (E-E e ) 

71=1 



(7) 



where a\ = dp uj (E e )/dE = X/N. While in the canonical ensemble there could 
exist three intersections points (e, c\ and c-i) because of the constancy of the 
inverse temperature, it is possible to ensure the existence of only one inter- 
section point E e by appropriately choosing the inverse temperature dependence 
j3 u (E). If the system size N is sufficiently large, the energy distribution function 
p (E) adopts a bell shape, which is approximately described with the Gaussian 
distribution: 

p(£)~Aexp --L(£-£; e ) 2 , (8) 

where oe = AS. The expectation values (E) and (/3 U ) provide a direct estima- 
tion of the energy E e and inverse temperature /3 e = (3 (E e ) at the intersection 
point illustrated in FIGHJ 



E e = (E) , f3 e = (fi u ) 



(9) 



The procedure previously described is equivalent to the one employed by 
Gerling and Hiiller in the framework of the dynamical ensemble [12] : 



uj d (E) — A (Et — E) 1 



(10) 
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which accounts for the thermal coupling of the system with a bath exhibiting a 
constant heat capacity, e.g., an ideal gas, whose inverse temperature obeys the 
following dependence on the system energy E: 

p u (E) = Bj (E T - E) . (11) 

Our proposal considers some improvements to the Gerling and Huller's method- 
ology. In fact, the energy dependence (fTTj) is less convenient for calculations than 
the linear ansatz ([3]), which is a feature particularly useful to perform the analy- 
sis of finite size effects (see subsection l2.2l belowV Moreover, one can also obtain 
the value of the heat capacity C (E), or more exactly, the so-called curvature 
curve k (E): 

K (E) = fN/C = -N^± (12) 

at the intersection point E ei n e = k (E e ), through the energy dispersion AE in 
analogous way of MC calculation by using the canonical ensemble (U} : 

Ke = (13 ) 
(AE) /N y ' 

This last equation was obtained from Eq. (fT2"| by rewriting the first relation of 
Eq.flU). 

Since the energy and the inverse temperature dispersions in Eq.Q are con- 
trolled by the coupling constant A, it is desirable to reduce them as low as 
possible. One can verify that the energy dispersion AE decreases with the in- 
crease of the coupling constant A. However, the value of this parameter should 
not be excessively large. Its increment leads to the increase of the inverse tem- 
perature dispersion A/3 U , which affects the precision of the inverse temperature 
f3 e of the system indirectly derived from the expectation value (f3 u ). A thermo- 
dynamic criterion to provide an optimal value for the coupling constant A can 
be obtained by minimizing the total dispersion A^: 

A| = ^+iV(AA,) 2 = i±f, (W) 
iV A + K e 

which talks about the precision in the determination of the intersection point 
(E e ,(3 e ). This analysis leads to the following result: 



Aa = Aa (« e ) = yl + K l - K e, minA^ = 2AA (15) 
2.2 Finite size effects 

The energy distribution function corresponding to the thermal coupling of the 
system with an environment can be expressed by the following equation: 

p{E)dE = w(E)n(E)dE. (16) 
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Here, ui (E) is a probabilistic weight that characterizes such a thermodynamic 
influence, which is related to the environment inverse temperature /3 U (E) as 
follows: 

(it) 

A direct integration allows us to verify that the probabilistic weight associated 
with the linear ansatz ([3]) is simply the Gaussian ensemble [15) : 



WG (E) 



1 



Zx(f3 e 



■ exp 



-PeE 



2iV^ 



E P 



introduced by Hetherington [TB] , which approaches in the limit A 
microcanonical ensemble: 

ujm = ^S(E- E e ) . 



(18) 
+oo to the 

(19) 



As already explained, the estimation of the microcanonical caloric curve 
(i (E) , as well as the fluctuation relation ((2|) , are based on the consideration of a 
Gaussian shape for the energy distribution function p (E) . Such an approxima- 
tion naturally arises as the asymptotic distribution as long as the system size 
N be sufficiently large. If the system size N is not so large, small deviations 
from the Gaussian profile © are naturally expected. Fortunately, the particular 
mathematical form of the Gaussian ensemble (fT8]) possibilities to consider some 
simple corrections formulae to deal with the existence of these finite size effects 
and to improve the precision of the present methodology. 

Let us denote the energy distribution function in terms of the energy per 
particle e = E/N as follows: 



p (e) = i- exp [-N</> (e)] 



(20) 



where the function (j> (e) = /3 e e+^X (e — e e ) 2 — s (e), with s (e) being the entropy 



per particle, and j3 e 



ds (e e ) 



(21) 



the system inverse temperature at the stationary point e e . Let us now develop 
a power-expansion: 



<t> (e e + x) — 4> (e e ) + - (A + Ke) x 2 + } a n x 



n— 3 



(22) 



with K e = —d 2 s (e e ) /de 2 being the curvature. The Gaussian approximation is 
developed by dismissing those terms with n > 2. Here, the energy deviation x 
obeys a Gaussian distribution: 



,(0) 



(a;) dx 



1 



2vrcr 



exp 



x 
'2^ 



dx 



(23) 



G 



with standard deviation a: 



" 2 " F(aW (24) 

and the partition function Z\ can be approximated as: 

Z x = e- N ^V2^(T. (25) 



Note that Eq. ([24|) is fully equivalent to Eq. (fT3[) since the standard deviation 
a = AE/N. Eq. ([25|) can be rewritten as follows: 

P x ~ /3 e £ e - S (E e ) + X - log (2tr7 2 ) , (26) 

where P\ = — log can be referred to as the Plank thermodynamic potential 
corresponding to the Gaussian ensemble fj 18[) . By performing the thermody- 
namic limit N — > oo, this thermodynamic function can be related to the known 
Legendre transformation with the microcanonical entropy: 

p 

PX = Jim -rr = fieSe ~ 9 (e e ) . (27) 
iv — »oo i\ 



Clearly, the Gaussian ensemble (|18|) provides a suitable extension of the Gibbs 
canonical ensemble ([TJ that is able to deal with the existence of macrostates 
with negative heat capacities C < 0. In fact, it is a particular example of the 
so-called generalized canonical ensembles that preserve some of its more relevant 
features [nJ[T7l[T8]. 

The first correction of the Gaussian approximation (1231) is obtained by dis- 
missing terms with n > 3 in the power expansion 



Ncj) (x) = N(/> + ^x 2 + £x 3 + O (x 4 ) , (28) 

where: 

« - -X^- 

It is convenient to introduce the dimensionless variable 8 — xj a. By considering 
the size dependencies £ ~ N and er ~ l/y/N, it is possible to verify that the 
cubic term £x 3 = ^a 3 9 3 decreases as l/y/~N with the increase of the system size 
N. Thus, one arrives at the distribution function: 

P w (e) aie ~ ^=e~^ (i - £o- 3 e 3 ) dd. (30) 

V 27T 

This last result leads to the following expectation values: 

(6) = -3£o- 3 , (6 2 ) = l, (9 3 ) = -15t<j 3 , (31) 
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which allow to express the expectation value of the energy deviation (x) as well 
as its second and third order dispersions, (fa 2 ) and (fa 3 ), as follows: 

(x) = - 3 ^ 4 + o(^),(fa 2 )=a 2 + o(^, (32) 

(fa 3 ) = -6^ + 0^. (33) 

The previous results can be combined with the linear ansatz to obtain the 
following first-order correction of Eq. © : 

* = <^- 2^ (^) 

^ = { ^- X 2Njsm^- (35) 

The second-order correction of the Gaussian constributions is carried out by 
dismissing terms with n > 4 in power expansion (1221) : 

N<P (x) = ^x 2 + fa 3 + 6x 4 + O (x 5 ) , (36) 

where: 

These terms lead to the following correction of the distribution function: 



lfl2 



j4\/27r V 2 



whose third and fourth terms account for finite size effects of order 0(1/N). 
By denoting the auxiliary constants C± and C2 as follows: 

d = fr 3 , C 2 = £ 2 <7 4 , (39) 
direct calculations allow to obtain the normalization constant A: 

A = 1 - 3C 2 + y C 2 , (40) 
as well as the following expectation values: 

(9) = -3d, (e 2 ) = 1 - 12C 2 + 45C 2 , (41) 

(6» 3 ) = -15Ci, (<9 4 ) = 3-102C 2 + 465C* 2 . (42) 

While the expressions (|3"2")l and (j3"3")l for the expectation values (x) and (fa 3 ) 
remain invariable under the second-order approximation, the second and fourth 
order dispersions (fa 2 ) and (fa 4 ) exhibit the following corrections: 

(fa 2 ) = a 2 (1 - 12C 2 + 36C 2 ) + O ( -L ) , (43) 



(fa 4 ) = (3 - 102C 2 + 339C 2 ) a 4 + O ( -L ) • (44) 



JV :i 
TV 4 
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By introducing the cumulants t\ and £2: 

the auxiliary constants G\ and C2 can be expressed as follows: 

C*? = — ei , C* 2 = — £2 + — £1 (46) 
1 36 ' 10 360 v ' 

Thus, the main work equations © and (| 13[) can be expressed in this second- 
order approximation as follows: 

(AE) /N y ' 

where ip\ is a second-order term defined by the cumulants t\ and £2 as: 

^jje 2 + He, (50) 

It is easy to verify that the surviving finite size effects of the above correction for 
the caloric curve (3 in terms of the energy per particle e are of order O (l/N 3 ), 
while the ones corresponding to the curvature curve k versus e are of order 
0(1/N 2 ). 

As a by-product of the previous analysis, one can obtain the third and fourth 
derivatives of the entropy per particle s (e) as follows: 



4 d\s(e e ) N 3 



where tp2 is another second-order term defined by cumulants £1 and £2 as: 

12 41 

"02 = y£2 + y^£l- (53) 

The underlying finite size effects of the dependence (3 versus e estimated by 
Eq. (|5"Tj) are of order O (l/N 2 ), while the ones corresponding to the dependence 
£4 versus e estimated by Eq. (|52"f are of order O (l/N). 
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2.3 Implementation 

As already commented in the introductory section, a general way to extend a 
particular canonical MC algorithm with transition probability W (X, — > Xf; (3b) 
by using the present methodology is achieved by replacing the inverse temper- 
ature (3b of the canonical ensemble (|TJ) with a variable inverse temperature, 
/3b — > Pu (E). The fulfilment of the detailed balance condition: 

Pu (X z )W(Xi-> X 3 ;0t)=p u (X, ) W (X, ^XnPi) (54) 

demands to use a certain value of the environment inverse temperature fts 
for both the direct and reverse process defined from the condition: 

^ff- = exp (#,6Ek,) , (55) 
Pu (Xf) 

which is hereafter referred to as the transition inverse temperature. Here, 
p u (X) = u[E (X)] represents the distribution function associated with the en- 
vironment inverse temperature f3 u (E), and 5Eif = Ej — Ei, the energy varying 
of the system during the transition, where Ei = E (Xi) and Ef — E (Xf). The 
case of the linear ansatz ([3]) is special, since its transition inverse temperature 
fil, is simply given by: 

PL = \ (ft, + Pi), (56) 

where /3£, and /?£ are the bath inverse temperatures at the initial and the final 
configurations respectively, = j3 u (Ei) and /?£ = j3 u (Ef). 

The direct applicability of the above result is restricted due to the final 
configuration Xf must be previously known in order to obtain the exact value 
of the transition inverse temperature . While such a requirement can be 
always satisfied in a local MC such as Metropolis importance sample [TH [5U] > 
the final configuration Xf is a priori unknown in non-local MC methods such 
as clusters algorithms [21j[2l[23[2l[25l[2i[27l[28l[29]. For these cases, one is 
forced to employ an approximated value of the transition inverse temperature 
/3* , e.g.: the inverse temperature at the initial configuration, = j3 u (Ei). 
Although the resulting MC algorithm does not obey the detailed balance, we 
have shown that the deviation of the asymptotic distribution function p u (X) 
from the exact distribution function p u (X) can be disregarded for N sufficiently 
large. This is the reason why this method is particularly useful to overcome slow 
sampling problems in large scale MC simulations. 

The use of an approximated value for the transition inverse temperature /3* 
is not longer appropriated when one is also interested in the study of systems 
with size N relatively small. Such an approximation introduces uncontrollable 
finite size effects that cannot be dealt by using work equations (|47H49[) . In this 
case, it is necessary to fulfil the detailed balance condition (|54|) to obtain the 
Gaussian profile (fT8| as asymptotic distribution function. The most general way 
to achieve this aim is to introduce a posteriori acceptance probability Wi^ / : 

w^ f = min 1 1, exp {-pJEif) j (57) 



10 



to accept or reject the final configuration Xf. Here, the terms: 

Wi^ f = W [X t -> Xf-pl) , Wf-H = W [X f -> X 4 ; /3£] (58) 

represent the transition probabilities of the direct and the reverse process, re- 
spectively, which should be calculated for the given canonical clusters MC algo- 
rithm. 



2.4 Iterative schemes 

Given a certain dependence of the environment inverse temperature (3^ (E), 
one can obtain from a MC run a punctual estimation of the system inverse 
temperature ft, the curvature m, as well as the third and fourth derivatives of 
the entropy per particle (f and at the i-th intersection point £j. These values 
can be considered to provide the next dependence /3w (£*). The linear ansatz 
([3]) can be rewritten in terms of the energy per particle as follows: 

ft, = ft +Ai (e-ej), (59) 

where e* and ft* are roughly estimations of the correct values and ft, which 
are employed here as seed parameters. The coupling constant Xi is provided by 
its optimal dependence (|T5|) : 

A 4 = A A «) (60) 

by using an estimation k* of the curvature at the intersection point £i . The seed 
values (e*, /3*, k*) can be obtained from the previous estimated values {si, 
by using the power expansions: 

e* +1 = Si + Ae (61) 

= ft - Ki Ae + ^CfAe 2 , (62) 

k* +1 =Ki- Cf Ae, (63) 

where Ae is the energy step. Here, it is convenient to consider the power- 
expansion up to the third derivative of the entropy per particle Cf , since the 
calculation of the fourth derivative £f can be only obtained with a sufficient 
precision after performing a very large MC run. 



2.5 Efficiency 

The efficiency of MC methods is commonly characterized in terms of the so- 
called decorrelation time t, that is, the minimum number of MC steps needed to 
generate effectively independent, identically distributed samples in the Markov 
chain [8]. Its calculation in this approach can be performed by using the ex- 
pression: 

t = lim tm = Hm — — - — ^ £m ^ (64) 
A/->oo a/^oo var (El) 
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where var(eAi) — (e 2 M ) — (sm) 2 is the variance of Em, which is defined as the 
arithmetic mean of the energy per particle e over M samples (consecutive MC 
steps): 

1 M 

However, the decorrelation time r only provides a partial view of the efficiency 
in the case of the extended canonical MC methods discussed in this work. In 
general, the efficiency is more appropriately characterized by the number of MC 
steps S needed to achieve the convergence of a given run. The question is that 
the number of MC steps S needed to achieve a convergence of the expectation 
value (x) of a given observable x also depends on its thermal dispersion Ax. For 
example, to obtain an estimation of the expectation value (x) with a statistical 
error e x < a, the number of MC steps S should obey the following inequality: 

S >I<^ (66) 

While the fluctuating behavior of a given observable x is an intrinsic system 
feature in canonical MC methods, this is not longer valid in the present frame- 
work. In this case, the fluctuating behavior crucially depends on the nature of 
the external influence acting on the system, e.g.: the thermal dispersion of the 
system energy AE depends on the coupling constant A in Eq.Q. In particu- 
lar, the number of MC steps S needed to obtain a point of the caloric curve 
(e, 0) with a precision ^/e| + e| < a should be evaluated in terms of the total 
dispersion introduced in Eq. ([T3|) as follows: 

s = # < 67 > 

We refer to the quantity rj = rA^ as the efficiency factor. Clearly, an extended 
canonical MC method is more efficient as smaller is its efficiency factor r\. 



3 Applications 

3.1 Potts model and its extended canonical MC algorithms 

For convenience, let us reconsider again the model system studied in our previous 
paper pQ: the g-state Potts model [S]: 

^EO-W). (68) 

m 

defined on the square lattice Lx L with periodic boundary conditions. The sum 
{ij} is over nearest neighbor sites, with <n = 1, 2, ... q being the spin variable 
on the i-th site. This model undergoes a continuous phase transition when 
q = 2 — 4, which turns discontinuous for q > 4. 
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The direct way to implement an extended canonical MC simulation of this 
model is by using the Metropolis importance sample [191 120) , whose transition 
probability is given by: 

W (A, -> Xf,pl) = min {1, exp {-^JE if ) } . (69) 

Besides this last local algorithm, the MC study of the model (|68p can be also 
carried out by using non-local MC methods such as the known Swendsen-Wang 
[2"T1 |2"2"] or Wolff's [53] clusters algorithms. As discussed elsewhere [5], such 
clusters algorithms are based on the consideration of the Fortuin-Kasteleyn 
theorem [SO ] 131 ] : 

z= e-f*' = £ /(I-P)^ 6 ^, (70) 

spins bonds 

which possibilities a mapping of this model system to a random clusters model 
of percolation. Here, p = 1 — e~^ s is the acceptance probability of bonds, N c 
is the number of clusters, b is the number of bonds, and Nb is the total number 
of possible bonds. 

In order to implement the extended versions of these clusters algorithms, let 
us denote by pt = 1 — e~@" the acceptance probability of bonds by starting from 
the initial configuration Xj. The transition probability of the direct process 
Wi-yf can be expressed as follows: 

W^ f = P \«{l- Pi ) b » +b «. (71) 

where b a is the number of inspected bonds which have been accepted, while 
bp + bd is the number of inspected bonds which have been rejected. At this 
point, it is also important to identify the number of rejected bonds bd which 
have been destroyed in the final configuration A/, as well as the number b c 
of created bonds. Note that these bonds are responsible of the energy varying 
SEif after the transition, SEif = bd — b c . In the reverse process, the accepted 
bonds b a of the direct process are also accepted with probability P f — 1 — e - ^, 
while the created bonds b c as well as the rejected bonds b p that have not been 
destroyed in the final configuration Xf are now rejected with probability 1 — pj. 
Thus, the transition probability of the reverse process can be expressed 

as follows: 

W f ^=p b f(l- Pf p +b *. (72) 

Within the canonical ensemble, where pi — Pf = 1 — e - ' 3 - 5 , it is easy to see that 
transition probability obeys the detailed balance condition: 

Wf ^ 1 -expipBSEij), (73) 



For a general case with A > 0, one should introduce the a posteriori acceptance 
probability ([57| in order to fulfil the detailed balance, which can be expressed 
as follows: 

«*-, = {!, exp (ft,)}, (74) 
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Figure 2: Microcanonical caloric and curvature curves associated with the 2D 
seven-state Potts model obtained from the application of the extended canonical 
MC algorithms (MIS: Metropolis importance sample, SW: Swendsen-Wang, 
and WF: Wolff's clusters algorithm). Such results can be compared with the 
ones obtained from a direct numerical differentiation of the entropy per particle 
s = log /N obtained by using the Wang-Landau sampling method. 

where the argument 9if depends on the integer numbers {b a ,b p ,b c ,bd)'- 

Bit = b a log (j-J - (b d - b c ) (2b p + b c + b d ) . (75) 
3.2 Results and discussions 

The results derived from the extended versions of the Metropolis importance 
sample, as well as the Swendsen-Wang and Wolff's clusters algorithms are shown 
in FIG[2]for the particular case of the 2D seven-state Potts model with L = 25. 
Each point of these dependencies have been obtained from MC runs with 10 6 
steps. For comparison purposes, we have also carried out the calculation of the 
caloric /3 (e) and the curvature k (e) curves by performing a direct numerical 
differentiation of the entropy per particles s(e) = log Q(E)/N obtained from 
the Wang-Landau sampling method |35j . which is shown in the inset panel 
(N = L 2 ). Although there exist a good agreement among all these MC results, 
the ones obtained from the Wang-Landau method seem to be less significant, 
overall, the results corresponding to the curvature curve k(e). 
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Figure 3: Microcanonical caloric curves j3 (e) showing the different corrections of 
finite size effects corresponding to the Swendsen-Wang clusters algorithm. We 
also show here the expectation value of the acceptance probability p = (m,-_»/) 
versus energy per particle s corresponding to the extended Swendsen-Wang 
algorithm obeying detailed balance. 

The previous results constitute a clear illustration that the finite size correc- 
tions formulae introduced in subsection 12.21 provide a significant improvement 
in the precision of this kind of MC calculations. In this particular example 
with L = 25, the contribution of first-order correction of finite size effects in the 
caloric curve j3 (e) have a typical order of 8\ ~ 1CP 3 , while the one corresponding 
to the second-order correction are of order <5 2 ~ 10~ 5 . Nevertheless, the most 
significative correction of finite size effects in non-local canonical MC methods 
comes from the fulfilment of the detailed balance condition (|54|) obtained after 
the introduction of a posteriori acceptance probability (l57|) . whose correction 
has a typical order of Sdb ~ 1CP 2 . This fact can be appreciated in FIG E] for the 
case of the Swendsen-Wang clusters algorithm. 

Although the acceptance probability for clusters flipping is lower than 

the unity, its expectation value p = (w^f) is significantly high for most of 
energy region (see also in FIG|3]). Since the consideration of such an a posteriori 
acceptance probability Wi-^f corrects a finite size effect error 8p = 1/3*, — /3*J 
associated with the estimation of the transition inverse temperature (3^, one 
should observe a growth of the expectation value (w^f) as N increases. Such 
a behavior is indeed appreciated in FIG H] for the extended clusters algorithms. 
While all these dependencies have a similar qualitative behavior, the expectation 



15 



1.0 



0.9 



0.8 





2D 7-state Potts model 




° L=8 


^^^^^^ 




— °— L=16 






—a— L=32 






—v— L=64 


extended SW 
■ .... ■ 





0.7 




0.0 0.2 0.4 0.6, 0.8 1.0 1.2 



extended WF 



0.0 0.2 0.4 



0.6 0.8 



1.0 1.2 



Figure 4: Size dependence of the expectation value of the a posteriori the ac- 
ceptance probability p — (w^f) for the extended clusters algorithms. 



value of the acceptance probability (wi^j) of the extended Wolff's clusters 
algorithm in the paramagnetic region is larger than the one corresponding to the 
extended Swendsen-Wang method. This is a quite expected qualitative result. 
The acceptance probability u>i_>./ should correct in the first case the finite size 
error Sp — | — | related to the flipping of only one cluster, while the finite 
size error 8p related to the extended Swendsen-Wang method is larger due to 
this algorithm involves the flipping of all system clusters. At low energies or in 
the ferromagnetic region, these methods have practically the same performance, 
since the number of spins belonging to the cluster is comparable to the system 
size N (see a more detailed comparison in the inset panel of FIG|4]). 

The qualitative behavior of the dependence (w^f) versus e finds a simple 
explanation in terms of the explicit dependence of the acceptance probability 
uii^ff on the coupling constant A. While the acceptance probability = 1 

within the canonical ensemble where A = 0, this quantity undergoes a reduction 
with the increase of the coupling constant A. Moreover, the optimal value of 
coupling constant Aa = vT+k? — k employed in the present MC simulations 
increases with the reduction of the system curvature k. Indeed, the lowest 
values of (w^f) are observed in the energy region where the system curvature 
also exhibits its lower values, that is, the anomalous region with negative heat 
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Figure 5: Main panel: Size dependence of the decorrelation time r and the 
efficiency factor -q — rAf, for the extended canonical MC simulations with 
environment inverse temperature /^(e) = /3 e +\(s— e e ). Here, /3 e = /3 pt , e e = £2, 
A = Aa (^2) and K2 — k (£2), where f3 pt is the estimated inverse temperature that 
corresponds to the discontinuous PT, while £2 is the stationary solution derived 
from the thermal equilibrium condition ^(£2) = Ppt located within the region 
with negative heat capacities C < 0. Inset panel: Size dependence of the entropy 
per particle s(e), or more exactly, the quantity s* (e) — s(e) — {3 p t£+ const, which 
allows to show the convex intruder associated with the region with negative 
heat capacities. The energies corresponding to the two maxima (£i,£3) and 
the minimum £2 are the three stationary solutions derived from the thermal 
equilibrium condition /? (£1,2,3) = fi p t- 



capacities C < 0. 

The size dependencies of the decorrelation time r and the efficiency factor 
77 = rAy within the anomalous region with C < are shown in the main panel 
of FIGjSJ Due to computational limitations, the data were obtained from MC 
simulations with lattice sizes L = (8, 16, 32, 64, 128). Since the absolute values 
of the curvature k are close to zero |k| ~ 0, the total dispersion is approximately 
given by the constant value Af^ ~ 2, in accordance with Eq. (fl~5|) . This is the 
reason why the dependencies t(N) and rj(N) are almost displaced in a constant 
value along the vertical direction of the log- log graph shown in FIG(5j 

As expected, the extended version of the Metropolis importance sample ex- 
hibits the largest values of the decorrelation time r for the lattice sizes L studied 
in this work. Curiously, the size dependence of its decorrelation time t(N) shows 



17 



an abrupt transition from a power-law regime with exponent a% ~ 0.72 to other 
one with exponent ai ~ 0.14 close to N ~ 10 3 (L — 32). Despite of the extended 
Metropolis importance sample is a local MC method, the effective exponent <xi 
for the larger system sizes is comparable to the ones associated with the ex- 
tended clusters MC methods, asw — 0.18 (Swendsen-Wang) and awF — 0.15 
(Wolff). Note also that the efficiency of these extended clusters methods are still 
significant despite of the consideration of the a posteriori acceptance probability 
([74]) employed here to ensure the detailed balance condition ([54]). 

The above examples confirm that the efficiency achieved with the application 
of the present methodology is more significant than the improvement considered 
by the application of re-weighting techniques such as multicanonical method 
and its variant, whose typical values of the exponent a ranges from 2 to 2.5 in 
the case of Potts models [34]. As already evidenced in FIG[2] while the Wang- 
Landau sampling method provides a good estimation of the entropy per particle 
s(e), its underlying statistical errors are still appreciable in the curvature curve 
k (s) = — d 2 s (e) jde 1 regardless the simulation was extended until the modifying 
factor fulfills the condition /j < exp(l0 -10 ). Such an observation evidences 
that the results obtained from this last method are not sufficiently relaxed to 
provide a precise estimation of the curvature curve «(e). By considering the 
CPU time-cost needed to achieve the convergence, it is more convenient to 
perform a punctual estimation of the caloric /3 (e) and curvature n (e) curves 
with any extended canonical MC method instead of carrying out a numerical 
differentiation of the microcanonical entropy s(e) over an energy range obtained 
from a re-weighting technics such as the Wang-Landau sampling method. 

4 Conclusions 

In this work, we have shown that the methodology to extend canonical MC 
methods inspired on the consideration of the recently obtained fluctuation re- 
lation ([2]) can be improved to account for the existence of finite size effects 
and to fulfil of the detailed balance condition ([54]) . Remarkably, despite of the 
consideration of the a posteriori acceptance probability (|57[) reduces the effi- 
ciency of the clusters MC methods, it has been shown that the relaxation times 
needed to ensure the convergence are more significant than the ones achieved 
with re- weighting technics such as multicanonical methods and its variants. For 
the particular case of the seven-states Potts model, the consideration of any 
extended canonical MC algorithm enables a suppression of the super-critical 
slowing down associated with the occurrence of the temperature driven first- 
order phase transition of this model from an exponential growth to a very weak 
power-law dependence with exponent a = 0.14 — 0.18. 

There is still some open questions in regard to the potentialities of the present 
methodology. For example, although this method has been specially conceived 
to overcome the slow relaxations of canonical MC simulations near to a region 
with a first-order phase transition, in principle, there is no limitation that the 
same one can be also employed to improve canonical MC simulations near to a 
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critical point of a continuous phase transition. Moreover, a similar extension is 
possible to carry out for those MC methods based on the consideration of the 
Boltzmann-Gibbs distributions: 

dp BG (E, X) = -^Y) ° XP [ ~^ E + YX)] dEdX (76) 

to account for the existence of anomalous values in other response functions 
besides the heat capacity. The analysis of these questions deserve a special 
attention in future works. 
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